/*
 * Copyright (c) Tecplot, Inc.
 *
 * All rights reserved.  Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 *   - Redistributions of source code must retain the above copyright notice, this list of
 *     conditions and the following disclaimer.
 *   - Redistributions in binary form must reproduce the above copyright notice, this list
 *     of conditions and the following disclaimer in the documentation and/or other
 *     materials provided with the distribution.
 *   - Neither the name of the Tecplot, Inc., nor the names of its contributors may be
 *     used to endorse or promote products derived from this software without specific
 *     prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
 * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
 * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
 * THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF
 * THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
#include "ComsolConverter.h"
#include "Status.h"
#include "TecFile.h"
#include <fstream>
#include <utility>
#include <algorithm>
#include <stack>

#ifdef _MSC_VER
#pragma warning (push)
#pragma warning (disable : 4512) // assignment operator could not be generated
#pragma warning (disable : 4127) // conditional expression is constant
#pragma warning (disable : 4701) // potentially uninitialized local variable 'result' used
#endif

#include "boost/algorithm/string/trim.hpp"
#include "boost/lexical_cast.hpp"

#ifdef _MSC_VER
#pragma warning (pop)
#endif

using namespace boost;
using std::string;
using std::vector;

// static member definition
ComsolConverter::ZoneMap_t    ComsolConverter::m_zoneTypeMap;
ComsolConverter::SectionMap_t ComsolConverter::m_stringToSectionMap;

/**
 */
ComsolConverter::Error::Error(const std::string& message):
    m_errorMessage(message)
{
}

/**
 */
ComsolConverter::Error::~Error() throw()
{
}

/**
 */
const char* ComsolConverter::Error::what() throw()
{
    return m_errorMessage.c_str();
}

/**
 */
void ComsolConverter::initZoneMap()
{
    if (m_zoneTypeMap.empty())
    {
        m_zoneTypeMap.insert(std::make_pair(ComsolZoneType_Points,      ZoneType_Ordered));
        m_zoneTypeMap.insert(std::make_pair(ComsolZoneType_Lines,       ZoneType_FELineSeg));
        m_zoneTypeMap.insert(std::make_pair(ComsolZoneType_Triangular,  ZoneType_FETriangle));
        m_zoneTypeMap.insert(std::make_pair(ComsolZoneType_Tetrahedral, ZoneType_FETetra));
        CHECK(m_zoneTypeMap.size() == END_ComsolZoneType_e);
    }
}

/**
 */
void ComsolConverter::initSectionMap()
{
    if (m_stringToSectionMap.empty())
    {
        m_stringToSectionMap.insert(std::make_pair(string("Times"),       ComsolSection_TimeParamEigenvalue));
        m_stringToSectionMap.insert(std::make_pair(string("Parameter"),   ComsolSection_TimeParamEigenvalue));
        m_stringToSectionMap.insert(std::make_pair(string("Eigenvalues"), ComsolSection_TimeParamEigenvalue));
        m_stringToSectionMap.insert(std::make_pair(string("Coordinates"), ComsolSection_Coordinates));
        m_stringToSectionMap.insert(std::make_pair(string("Elements"),    ComsolSection_Elements));
        m_stringToSectionMap.insert(std::make_pair(string("Data"),        ComsolSection_Data));
    }
}

/**
 */
ComsolConverter::ComsolConverter(const std::string& fileName):
    m_fileName(fileName)
{
    /*
     * Make sure we can open the file for read access.
     */
    std::ifstream testStream;
    openFileForRead(fileName, testStream);
    initZoneMap();
    initSectionMap();
}

/**
 */
ComsolConverter::~ComsolConverter()
{
}

/**
 */
void ComsolConverter::convertToTecplot(const std::string& destinationFile)
{
    std::ifstream inputFile;
    openFileForRead(m_fileName, inputFile);

    std::vector<ZoneInfo_s> zoneInfoList;
    scanFile(inputFile,
             zoneInfoList);
    if (zoneInfoList.empty())
    {
        throw Error("No data found in file.");
    }

    vector<string> variables;
    createVariableList(zoneInfoList,
                       variables);

    try
    {
        TecFile outputFile("Comsol data set",
                           variables,
                           destinationFile,
                           ".", // scratch dir
                           TecFile::Float);

        vector<ZoneInfo_s>::iterator zoneIter = zoneInfoList.begin();
        Status<Error> status("Reading zone data...",
                             zoneInfoList.size());
        while (zoneIter != zoneInfoList.end())
        {
            vector<ZoneData_s> zoneData(m_solutionTimes.size());
            readZoneData(*zoneIter,
                         inputFile,
                         zoneData);
            outputZoneData(*zoneIter,
                           zoneData,
                           variables,
                           outputFile);

            ++zoneIter;
            ++status;
        }

        outputFile.close();
    }
    catch (TecFile::Error& e)
    {
        // Translate exception
        throw Error(string("Error creating Tecplot file:\n") + e.what());
    }
}

/**
 */
void ComsolConverter::openFileForRead(const std::string&   fileName,
                                      std::ifstream& stream) const
{
    // Open in binary mode to prevent issues with Unix line endings and seekg() on Windows
    stream.open(fileName.c_str(), std::ios::in | std::ios::binary);
    if (stream.fail())
        throw Error("Unable to open file for reading.");
}

/**
 */
void ComsolConverter::scanFile(std::ifstream&           inputFile,
                               std::vector<ZoneInfo_s>& zoneInfoList)
{
    REQUIRE(inputFile.good());

    inputFile.seekg(0, std::ios::end);
    Status<Error> status("Scanning file for data...",
                         inputFile.tellg());
    inputFile.seekg(0, std::ios::beg);

    ComsolSection_e curSection   = ComsolSection_Invalid;
    LgIndex_t       curItemCount = 0; // Keeps track of either the number of nodes or elements as we scan
    ZoneInfo_s      curZoneInfo;
    string          line;

    std::getline(inputFile, line);
    while (!inputFile.eof() &&
           inputFile.good())
    {
        trim(line);

        if (lineIsHeader(line))
        {
            ComsolSection_e newSection = determineHeaderType(line);
            switch (newSection)
            {
                case ComsolSection_TimeParamEigenvalue:
                {
                    // Solution times are on the next line
                    std::getline(inputFile, line);
                    readSolutionTimes(line, m_solutionTimes);
                } break;

                case ComsolSection_Coordinates:
                {
                    /*
                     * We've come across a new zone. Push the current zone info onto the infoList
                     * (if it's valid), and reset curZoneInfo struct.
                     */
                    if (curZoneInfo.zoneType != ComsolZoneType_Invalid)
                    {
                        if (curSection == ComsolSection_Elements)
                            curZoneInfo.numElements = curItemCount;

                        // Assign strandID for transient data
                        if (!m_solutionTimes.empty())
                            curZoneInfo.strandID = static_cast<Strand_t>(zoneInfoList.size() + 1);

                        zoneInfoList.push_back(curZoneInfo);
                    }

                    curZoneInfo.reset();
                    curItemCount = 0;

                    readCoordinateHeader(line, inputFile, curZoneInfo.zoneType);

                    curZoneInfo.nodeDataPosition = inputFile.tellg();
                    curSection = ComsolSection_Coordinates;
                } break;

                case ComsolSection_Elements:
                {
                    /*
                     * We've come to the end of the node values.
                     */
                    curZoneInfo.numNodes = curItemCount;
                    curItemCount = 0;

                    curZoneInfo.elementDataPosition = inputFile.tellg();
                    curSection = ComsolSection_Elements;
                } break;

                case ComsolSection_Data:
                {
                    /*
                     * We've come to a new data section. The data section(s) contain
                     * the solution data for the zone.
                     */
                    curZoneInfo.numElements = curItemCount;

                    std::string varName;
                    readDataHeader(line,
                                   varName);
                    curZoneInfo.varNames.push_back(varName);
                    curZoneInfo.solutionDataPositions.push_back(inputFile.tellg());

                    curSection = ComsolSection_Data;
                } break;

                default:
                    CHECK(FALSE);
            }
        }
        else if (!line.empty() &&
                 (curSection == ComsolSection_Coordinates ||
                  curSection == ComsolSection_Elements))
        {
            ++curItemCount;
        }

        std::getline(inputFile, line);
        status.update(inputFile.tellg());
    }

    if (!inputFile.eof() &&
        inputFile.fail())
    {
        throw Error("File I/O error.");
    }
    else
    {
        inputFile.clear();
        inputFile.seekg(0, std::ios::beg);
    }

    // Don't forget to push_back the last zone
    if (curZoneInfo.zoneType != ComsolZoneType_Invalid)
    {
        if (curSection == ComsolSection_Elements)
            curZoneInfo.numElements = curItemCount;

        // Assign strandID for transient data
        if (!m_solutionTimes.empty())
            curZoneInfo.strandID = static_cast<Strand_t>(zoneInfoList.size() + 1);

        zoneInfoList.push_back(curZoneInfo);
    }

    // Need at least one solution time
    if (m_solutionTimes.empty())
        m_solutionTimes.push_back(0.0);
}

/**
 */
ComsolConverter::ComsolSection_e ComsolConverter::determineHeaderType(const std::string& headerLine) const
{
    REQUIRE(lineIsHeader(headerLine));

    // Extract the keyword from the line
    string::size_type keywordStart = headerLine.find_first_not_of('%'); // Skip past the %
    keywordStart = headerLine.find_first_not_of(' ', keywordStart + 1); // Skip past intermediate spaces
    string::size_type keywordEnd   = headerLine.find_first_of(' ', keywordStart + 1);
    string keyword = headerLine.substr(keywordStart,
                                       keywordEnd - keywordStart);

    SectionMap_t::iterator sectionIter = m_stringToSectionMap.find(keyword);
    if (sectionIter == m_stringToSectionMap.end())
    {
        std::ostringstream errMessage;
        errMessage << "Unknown header found: " << keyword;
        throw Error(errMessage.str());
    }

    ENSURE(VALID_ENUM(sectionIter->second, ComsolSection_e));
    return sectionIter->second;
}

/**
 */
void ComsolConverter::readElementHeader(const std::string& headerLine,
                                        ComsolZoneType_e&  zoneType) const
{
    REQUIRE(lineIsHeader(headerLine));

    if (headerLine.find("triangular") != string::npos)
        zoneType = ComsolZoneType_Triangular;
    else if (headerLine.find("lines") != string::npos)
        zoneType = ComsolZoneType_Lines;
    else if (headerLine.find("points") != string::npos)
        zoneType = ComsolZoneType_Points;
    else if (headerLine.find("tetrahedral") != string::npos)
        zoneType = ComsolZoneType_Tetrahedral;
    else
    {
        std::string errMessage = "Unsupported element type. Error occurred on this line:\n";
        errMessage += headerLine;
        throw Error(errMessage);
    }
}

/**
 * Extract the variable name from the header line
 */
std::string ComsolConverter::readVariableNameFromHeader(std::string const& headerLine) const
{
    // The variable name is between parentheses
    std::string::const_iterator varNameStart = std::find(headerLine.begin(), headerLine.end(), '(');

    if (varNameStart == headerLine.end())
        throwMissingVariableNameError(headerLine);

    // Find the matching end parentheses. The name may contain a nested set of
    // parentheses, e.g. (sin(u)), so we use a stack
    std::stack<char> parenStack;
    parenStack.push(*varNameStart);
    ++varNameStart;
    std::string::const_iterator varNameEnd = varNameStart;
    while (!parenStack.empty())
    {
        ++varNameEnd;
        if (*varNameEnd == '(')
            parenStack.push(*varNameEnd);
        if (*varNameEnd == ')')
            parenStack.pop();
    }

    return std::string(varNameStart, varNameEnd);
}

void ComsolConverter::readDataHeader(const std::string& headerLine,
                                     std::string&       varName) const
{
    REQUIRE(determineHeaderType(headerLine) == ComsolSection_Data);

    varName = readVariableNameFromHeader(headerLine);
    trim(varName);

    if (varName.empty())
        throwMissingVariableNameError(headerLine);
}

/**
 */
bool ComsolConverter::lineIsHeader(const std::string& line) const
{
    return !line.empty() &&
           line.at(0) == m_comsolSectionDelim;
}

/**
 */
void ComsolConverter::readZoneData(const ZoneInfo_s&        zoneInfo,
                                   std::ifstream&           inputFile,
                                   std::vector<ZoneData_s>& zoneData) const
{
    REQUIRE(inputFile.good());
    REQUIRE(!zoneData.empty());

    readCoordinateData(zoneInfo,
                       inputFile,
                       zoneData);

    readConnectivityData(zoneInfo,
                         inputFile,
                         zoneData.at(0).connectivity);
    // Connectivity for other zones is the same, so just copy
    std::vector<ZoneData_s>::iterator zoneIter = zoneData.begin();
    ++zoneIter;
    while (zoneIter != zoneData.end())
    {
        zoneIter->connectivity = zoneData.front().connectivity;
        ++zoneIter;
    }

    readSolutionData(zoneInfo,
                     inputFile,
                     zoneData);
}

/**
 */
void ComsolConverter::readCoordinateData(const ZoneInfo_s&        zoneInfo,
                                         std::ifstream&           inputFile,
                                         std::vector<ZoneData_s>& zoneData) const
{
    REQUIRE(zoneData.size() == m_solutionTimes.size());
    REQUIRE(inputFile.good());

    const bool zoneHasThirdDimension = zoneInfo.zoneType == ComsolZoneType_Tetrahedral;

    // Resize each array to the number of nodes
    std::vector<ZoneData_s>::iterator zoneIter = zoneData.begin();
    while (zoneIter != zoneData.end())
    {
        zoneIter->xVals.reserve(zoneInfo.numNodes);
        zoneIter->yVals.reserve(zoneInfo.numNodes);
        if (zoneHasThirdDimension)
            zoneIter->zVals.reserve(zoneInfo.numNodes);
        ++zoneIter;
    }

    inputFile.seekg(zoneInfo.nodeDataPosition);
    CHECK(inputFile.tellg() == zoneInfo.nodeDataPosition);

    // Each line contains coordinates for a single node across all time steps
    for (LgIndex_t curNode = 0; curNode < zoneInfo.numNodes; ++curNode)
    {
        string line;
        std::getline(inputFile, line);
        std::istringstream lineBuffer;
        lineBuffer.str(line);

        std::vector<float> values;
        std::copy(std::istream_iterator<float>(lineBuffer), std::istream_iterator<float>(),
                  std::back_inserter(values));

        // The line contains either a single coordinate tuple (for static meshes), or a
        // coordinate tuple for each time step
        size_t numCoords = zoneHasThirdDimension ? 3 : 2;
        bool hasMultipleTimes = (m_solutionTimes.size() > 1) &&
                                (values.size() == m_solutionTimes.size() * numCoords);

        std::vector<float>::iterator valueIter = values.begin();
        for (size_t curTime = 0; curTime < m_solutionTimes.size(); ++curTime)
        {
            zoneData.at(curTime).xVals.push_back(*(valueIter++));
            zoneData.at(curTime).yVals.push_back(*(valueIter++));
            if (zoneHasThirdDimension)
                zoneData.at(curTime).zVals.push_back(*(valueIter++));

            if (!hasMultipleTimes)
                valueIter = values.begin();
        }
    }
}

/**
 */
void ComsolConverter::readConnectivityData(const ZoneInfo_s&   zoneInfo,
                                           std::ifstream&      inputFile,
                                           ConnectivityList_t& connectivity) const
{
    REQUIRE(inputFile.good());

    // We can skip connectivity for Point zones
    if (zoneInfo.zoneType == ComsolZoneType_Points)
        return;

    // Determine how many values to read per element
    size_t nodesPerElement = 0;
    switch (zoneInfo.zoneType)
    {
        case ComsolZoneType_Points:      nodesPerElement = 1; break;
        case ComsolZoneType_Lines:       nodesPerElement = 2; break;
        case ComsolZoneType_Triangular:  nodesPerElement = 3; break;
        case ComsolZoneType_Tetrahedral: nodesPerElement = 4; break;
        default:
            CHECK(FALSE);
    }

    connectivity.reserve(zoneInfo.numElements);
    inputFile.seekg(zoneInfo.elementDataPosition);
    CHECK(inputFile.tellg() == zoneInfo.elementDataPosition);
    for (LgIndex_t curElement = 0; curElement < zoneInfo.numElements; ++curElement)
    {
        for (size_t curNode = 0; curNode < nodesPerElement; ++curNode)
        {
            LgIndex_t nodeIndex = 0;
            inputFile >> nodeIndex;

            if (!inputFile.good())
            {
                string errMessage("Invalid data value found: ");
                string badValue;
                inputFile.clear();
                inputFile >> badValue;
                throw Error(errMessage += badValue);
            }

            connectivity.push_back(nodeIndex);
        }
    }
}

/**
 */
void ComsolConverter::readSolutionData(const ZoneInfo_s&        zoneInfo,
                                       std::ifstream&           inputFile,
                                       std::vector<ZoneData_s>& zoneData) const
{
    REQUIRE(zoneData.size() == m_solutionTimes.size());
    REQUIRE(inputFile.good());

    // Resize each array to the number of nodes. This is similar to what we do in
    // readCoordinateData, but for a variable number of items
    std::vector<ZoneData_s>::iterator zoneIter = zoneData.begin();
    while (zoneIter != zoneData.end())
    {
        zoneIter->solutionData.resize(zoneInfo.solutionDataPositions.size());
        std::vector<SolutionData_s>::iterator varIter = zoneIter->solutionData.begin();
        while (varIter != zoneIter->solutionData.end())
        {
            varIter->values.reserve(zoneInfo.numNodes);
            ++varIter;
        }
        ++zoneIter;
    }

    // Loop over variables, then nodes, then timesteps, since that's how things are stored in the file
    for (size_t curVar = 0; curVar < zoneInfo.solutionDataPositions.size(); ++curVar)
    {
        inputFile.seekg(zoneInfo.solutionDataPositions.at(curVar));

        // Each line contains values for a single node across all time steps.
        for (LgIndex_t curNode = 0; curNode < zoneInfo.numNodes; ++curNode)
        {
            string line;
            std::getline(inputFile, line);
            std::istringstream lineBuffer;
            lineBuffer.str(line);

            for (size_t curTime = 0; curTime < m_solutionTimes.size(); ++curTime)
            {
                float value;
                if (!(lineBuffer >> value))
                {
                    string errMessage("Invalid data value found: ");
                    string badValue;
                    lineBuffer.clear();
                    lineBuffer >> badValue;
                    throw Error(errMessage += badValue);
                }

                // Store the value.
                zoneData.at(curTime).solutionData.at(curVar).values.push_back(value);

            } // End loop over timesteps
        } // End loop over nodes
    } // End loop over variables
}

/**
 */
void ComsolConverter::readSolutionTimes(const std::string&  line,
                                        std::vector<float>& solutionTimes)
{
    std::istringstream input;
    input.str(line);

    while (!input.eof() &&
           input.good())
    {
        float solutionTime;
        input >> solutionTime;

        if (input.good())
            solutionTimes.push_back(solutionTime);
    }

    if (!input.eof() &&
        input.fail())
    {
        throw Error("File I/O error.");
    }
}

/**
 */
void ComsolConverter::readCoordinateHeader(const std::string& headerLine,
                                           std::istream&      inputFile,
                                           ComsolZoneType_e&  zoneType) const
{
    REQUIRE(determineHeaderType(headerLine) == ComsolSection_Coordinates);

    // Supported Comsol zones are either of element type triangle or tetrahedral, and we
    // can tell this based on the number of spatial coordinate variables
    size_t numCoordinateVars = 0;

    // The header will either contain "(x,y,z)" or "(x,y)", corresponding to Tetrahedral
    // or Triangular data. Easiest thing to do is count the number of commas
    size_t numCommas = std::count(headerLine.begin(), headerLine.end(), ',');
    if (numCommas > 0)
    {
        numCoordinateVars = numCommas + 1;
    }
    else
    {
        // Sometimes the header is missing the specification for the spatial dimensions.
        // In that case, just read the next line, which contains the first set of
        // coordinates, and count the numbers of values that can be read
        std::istream::pos_type fileLocation = inputFile.tellg();
        std::string nextLine;
        getline(inputFile, nextLine);
        inputFile.seekg(fileLocation, std::ios::beg);

        std::istringstream lineStream(nextLine.c_str());
        numCoordinateVars = std::distance(std::istream_iterator<double>(lineStream), std::istream_iterator<double>());
    }

    switch (numCoordinateVars)
    {
        case 2: zoneType = ComsolZoneType_Triangular; break;
        case 3: zoneType = ComsolZoneType_Tetrahedral; break;
        default:
        {
            std::string errMessage = "Invalid %Coordinates header. The number of spatial dimensions is invalid.\n";
            errMessage += "Error occurred on this line:\n";
            errMessage += headerLine;
            throw Error(errMessage);
        }
    }
}

/**
 */
void ComsolConverter::createVariableList(const std::vector<ZoneInfo_s>& zoneInfoList,
                                         std::vector<std::string>&      variables) const
{
    // Loop over all variables, determining how many spatial variables we need and how many
    // solution variables we have.
    std::vector<std::string> tempVariables;
    size_t maxDimensions = 0;

    std::vector<ZoneInfo_s>::const_iterator zoneIter = zoneInfoList.begin();
    while (zoneIter != zoneInfoList.end())
    {
        size_t zoneDimensions = 0;
        if (zoneIter->zoneType == ComsolZoneType_Triangular)
            zoneDimensions = 2;
        else if (zoneIter->zoneType == ComsolZoneType_Tetrahedral)
            zoneDimensions = 3;
        else
            throw Error("Error creating variable list.");

        maxDimensions = std::max(maxDimensions, zoneDimensions);

        // Add any new variables to our total variable list
        std::vector<std::string>::const_iterator varIter = zoneIter->varNames.begin();
        while (varIter != zoneIter->varNames.end())
        {
            std::vector<std::string>::iterator totalVarIter = std::find(tempVariables.begin(),
                                                                        tempVariables.end(),
                                                                        *varIter);
            if (totalVarIter == tempVariables.end())
                tempVariables.push_back(*varIter);

            ++ varIter;
        }

        ++zoneIter;
    }

    // Construct the final variable list
    variables.clear();
    variables.push_back("X");
    variables.push_back("Y");
    if (maxDimensions == 3)
        variables.push_back("Z");
    variables.insert(variables.end(),
                     tempVariables.begin(),
                     tempVariables.end());
}

/**
 */
void ComsolConverter::outputZoneData(const ZoneInfo_s&               zoneInfo,
                                     const std::vector<ZoneData_s>&  zoneData,
                                     const std::vector<std::string>& variables,
                                     TecFile&                        outputFile) const
{
    INVARIANT(zoneData.size() == m_solutionTimes.size());

    for (size_t curZone = 0; curZone < zoneData.size(); ++curZone)
    {
        outputFile.zone(string(), // zone name
                        m_zoneTypeMap[zoneInfo.zoneType],
                        zoneInfo.numNodes,
                        zoneInfo.numElements,
                        TecFile::Block,
                        zoneInfo.strandID,
                        m_solutionTimes.at(curZone));

        // Output spatial variables
        size_t firstNonSpatialVarIndex = 2;
        outputFile.data(zoneData.at(curZone).xVals);
        outputFile.data(zoneData.at(curZone).yVals);
        if (!zoneData.at(curZone).zVals.empty())
        {
            firstNonSpatialVarIndex = 3;
            outputFile.data(zoneData.at(curZone).zVals);
        }

        // Output non-spatial variables
        for (size_t curVar = firstNonSpatialVarIndex; curVar < variables.size(); ++curVar)
        {
            vector<string>::const_iterator varIter = std::find(zoneInfo.varNames.begin(),
                                                               zoneInfo.varNames.end(),
                                                               variables.at(curVar));
            if (varIter != zoneInfo.varNames.end())
            {
                size_t varIndex = std::distance(zoneInfo.varNames.begin(), varIter);
                outputFile.data(zoneData.at(curZone).solutionData.at(varIndex).values);
            }
            else
            {
                // We don't have data for this variable for this zone. Output all zeroes.
                outputFile.data(std::vector<float>(zoneInfo.numNodes, 0));
            }
        }

        // Output connectivity
        outputFile.connectivity(zoneData.at(curZone).connectivity);
    }
}

void ComsolConverter::throwMissingVariableNameError(std::string const& headerLine) const
{
    std::string errMessage = "Missing variable name in %Data header. Error occurred on this line:\n";
    errMessage += headerLine;
    throw Error(errMessage);
}